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ABSTRACT 

We  consider  plane-wave  motion  at  normal  inci- 
dence in  a horizontally  layered  system.  The  sys- 
tem is  assumed  lossless,  and  only  the  compres- 
sional  waves  are  treated.  A procedure  is  introduced 
for  determining  the  reflection  coefficients  of  the 
layered  system  when  the  observed  seismic  data 
may  contain  random  noise.  No  deconvolution  of  the 
measured  seismic  data  is  required  by  the  proce- 
dure when  the  input  is  a narrow  wavelet. 

1.  INTRODUCTION 

In  recent  years  much  attention  has  been  given  to 
the  problem  of  determining  reflection  coefficients 
for  a layered  media  from  the  observed  seismic^ 
data  [1-4].  In  line  with  the  customary  assumptions 
and  restrictions,  we  also  limit  our  attention  to  a 
horizontally  stratified  nonabsorptive  earth  with 
vertically  traveling  plane  compressional  waves. 
Such  a system  is  completely  described  by  a set  of 
reflection  coefficients  and  travel  times  within  lay- 
ers. 

A fundamental  procedure  described  in  detail  in 
the  above  references  for  deriving  values  of  the  re- 
flection coefficients  can  be  summarized  by  the  fol- 
lowing assumptions  and  steps. 

Standard  Assumptions: 

(Al)  The  input  wavelet  is  assumed  known. 

(A2)  The  data  is  assumed  noise  free. 

(A3)  The  layered  system  is  assumed  to  have  uni- 
form travel  times  between  layers  where  a number 
of  the  layers  are  hypothetical,  i.  e. , they  may  not 
correspond  to  an  actual  interface  of  the  substruc- 
ture and  are  associated  with  zero  reflection  and 
unity  transmission  coefficients. 

Standard  Steps; 

(SI)  The  observed  seismic  data  is  deconvolved  us- 
ing the  input  waveform  to  produce  the  system  re- 
sponse to  a unit  spike  input. 

tS2)  The  number  of  layers  is  chosen  high  enough  to 
result  in  travel  times  short  compared  with  the  in- 
verse of  the  bandwidth  of  the  observed  seismic 
data. 

(53)  The  deconvolved  data  is  sampled  with  samp- 
ling Interval  equal  to  the  chosen  one-way  travel 
time  between  lavers.* 

(54)  The  system  structure  is  used  to  arrive  at  a 
set  of  normal  equations  (linear  simultaneous  equa- 
tions) in  terms  of  reflection  coefficients  and  the 


discretized  and  deconvolved  observed  data. 

(S5)  The  normal  equations  of  the  preceding  step 
have  the  Toeplitz  structure  which  makes  it  possible 
to  utilize  the  very  efficient  Levinson  algorithm  to 
recursively  solve  for  the  reflection  coefficients. 

In  this  paper  the  method  of  solution  to  the  inverse 
problem  stated  above  is  fundamentally  modified  to 
cope  with  the  existence  of  the  noise  in  the  measure- 
ment data,  often  without  need  for  any  deconvolution. 
More  specifically,  although  again  a uniform  layered 
system  is  assumed,  the  choice  of  number  of  layers 
can  now  be  made  independent  of  the  sampling  rate 
requirement  of  the  data  (step  (S2)  above)  often  re- 
sulting in  the  need  for  far  fewer  layers.  No  decon- 
volution is  necessary  (step  (SI))  for  wavelets  of  dur- 
ation of  the  order  of  twice  the  layer  travel  times. 
The  exact  deconvolution  of  step(Sl)is  either  not 
possible  in  practice  or,  at  the  least,  will  further  ag- 
gravate the  harmful  effects  of  the  noise  in  the  ob- 
servation [S],  Furthermore,  the  deconvolution  is  a 
time  consuming  operation.  Finally,  the  procedure  is 
very  simple  to  derive  and  does  not  need  the  con- 
cepts of  z-transforms,  minimum  phase,  forward  and 
backward  polynomials,  spectral  factorization,  etc. 
The  results  reduce  to  the  existing  solution  of  the  in- 
verse problem  in  the  absence  of  noise  and  with  a 
spike  input  signal  (wavelet)  [l]. 

2.  STATEMENT  OF  THE  PROBLEM 

We  are  considering  a uniform  K layered  system 
and  normal  incident  compressional  waves.  Figure  1 
represents  such  a system  where  d,(t)  is  the  down- 
going wave  at  the  bottom  of  the  layer  and  u.(t)  is 
the  up-going  wave  at  the  top  of  the  layer.  The'  re- 
flection, downward  transmission  and  upward  trans- 
mission coefficients  associated  with  the  interface  at 
the  bottom  of  layer  are  denoted  rj,  t.  and  t',  re- 
spectively where  tj=l+r,,  tj'=l-rj.  Tile  one  way  trav- 
el time  between  layers  Is  denoted  by  r. 

The  input  to  the  system,  d^lt),  is  assumed  known 
(the  wavelet)  and  the  output  may  be  either  Uj(t)  (in 
the  marine  environment)  or  Ug(t).  The  measured 
seismic  data,  y(t),  consists  of  the  output  and  an  addi- 
tive noise  component  n(t).  The  source  of  this  noise 
may  be  the  instrument  measurement  noise,  the  un- 
certainty in  the  knowledge  of  the  input  wavelet  or 
response  to  unwanted  inputs  (ambient  noise).  It  is 
desired  to  process  y(t),ts0  and  derive  values  for 

the  reflection  coefficients  r.,  j = 1 Kir,  may  be 

assumed  known  in  cases  such  as  the  marine  envi- 
ronment). 


Of  course  the  deconvolution  may  be  performed  in 
discrete  time  using  the  same  sampling  interval. 


3.  STATE  EQUATIONS 

Using  the  notation  of  Flg.l,for  a general  layer  j 
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we  have  [6,  7]. 

Uj(t+T)  = + rjd_.(t)  (1) 

d..i(t+T)  = -rjUj^j(t)  +tjd,(t)  . (2) 

These  equations  are  valid  for  j=  1 K-1.  The 

They  should  be  augmented  at  the  surface  with 

u,(t)  = t'ui(t)+r,d,(t)  (3) 

d,(t+T)  = -roU,(t)  + t,d,(t)  (4) 

and  at  the  basement  with 

u,(t+T)  = tiu,*j(t)  + r,d,(t)  = r,d,(t)  (5) 

+ tidjit)  = t,d,(t)  . (6) 

Equations  (3,4)  and  (5,6)  can  be  derived  from  ( 1)  and 


(2)  [letting  j = 0,  1 K]  by  noting  that  u^lt)  is  talcen 

at  the  bottom  of  layer  0 and  d,^j(t)  represents  the 
down-going  wave  leaving  the  last  interface  and  la 
not  reflected  by  any  other  interface:  hence  u,^^(t)  aO, 
These  equations,  called  causal  functional,  are  not 
difference  equations  since  t is  the  continuous  time 
variable  [6], 

Using  the  state  equations  given  above,  it  can  be 
shown [Appen.  B,  9]  that  the  function  d,^j  satisfies 
the  equation 

d,*i[t+KT]  + a.d,.,[t+(K-2)T]+-  • • +a,.id,.i[t-(K-2)T] 

■'■*’ori(d,,i[t-KT]  = ^tidsd)  • (7) 

Note  that  the  coefficient  of  the  highest  term  of  the 
left  hand  aide  is  unity  and  that  of  the  lowest  term  is 
r^r,.  The  precise  form  of  the  other  coefficients  is 
not  important.  Only  the  structural  form  of  (7)  will 
be  utilized  in  the  sequel.  In  this  equation,  the  un- 
knowns are  the  reflections  coefficients  which  are 
embedded  in  the  coefficients  a^, , , . , a,.^,  r^r,  and 
^:t,.  The  input  dg(t)  is  assumed  known.  Equation (7) 

is  the  starting  point  for  our  inverse  procedure: 
however,  it  is  in  terms  of  a signal  which,  in  general, 
is  not  measurable.  In  the  following  two  sections  we 
relate  d,^^(t)  to  measured  seismic  data.  We  do  this 
so  that  we  will  be  able  to  extract  the  reflection  co- 
efficients from  measured  data. 

4.  A GENERALIZED  ENERGY  TRANSFER 

(KUNETZ)  RE1.ATIC5N 

Consider  e to  be  a noQoaegative  continuous  or 
discrete  variable  with  dimension  of  time«  Equations 
(1)  and  (2)  where  j=0,  1,  • • . , K are  multiplied  by 
u,(t+-+€)  and  dj^Jt+T+c)  respectively  resulting  in 

Uj(t+T)Uj(t+-+e)  = tJu.^,(t)Uj,j(t+e)+rJdj(t)dj(t+e) 

+Uj^j(t+e)dj(l)]  (8) 

dj»i(t+T)d.,j(t+T+e)  = r^J„^(t)u,^^(t+e)+tJdJ(t)d,(t+e) 

- rjtj[uj.j(t)dj(t+e)+Uj„j(t+£)d,(t)] . (9) 

Multiplying  (8)  by  t./tj  and  adding  the  resulting  ex- 
pression to  (9)  yields 

d..i(t+r)dj„j(t4-+e)+(tj/tpuj(t+T)uj(t+»+c) 

» (t,/t5)dj(t)dj(t+£)+u,„^(t)u.^j(t+€) . (10) 

Let  us  define  the  following  correlation-type 
functions 

D;(£)  dj(t)d,(t+s)dt  (11) 


Uj(e)  ^ j’^“uj(t)upt+e)dt.  (12) 

Integrating  both  sides  of  Eq.(lO)  from  -•  to  +•,  and 
using  Eqs.  (1 1)  and  (12),  we  find  that 

°j*i<e)-U,*i(e)  = t,/tj[Dj(e)-Uj(£)]  (13) 

where  j = 0,  1,  2, , . , , K.  This  is  a generalization  of 
the  well-known  [l]  energy  transfer  (Kunetz)  relation. 
Note  that  in  ourderivation,inputdg(t)  is  not  assumed 
to  be  an  impulse  and  the  seismic  data  is  not  dis- 
cretized. 

Iterating  ( 13),  starting  with  and  ending  with 
j=K,  we  obtain 

D,*i(e)  = -^[D^(c)-Uj(€)]  (14) 

where  t can  take  on  the  values  of  0,  1, ... , or  K, 
i A |i 

and  = In  the  marine  case  this  rela- 
0 0 

tionship  is  used  with  X = l.  In  the  non-marine  case 
it  is  used  with  X=0. 

5,  APPLICATION  TO  MARINE  ENVIRONMENT 

In  this  section  we  will  direct  our  attention  to  the 
marine  case  and  shall  express  0,,.i(e)  in  terms  of 
measured  signals.  To  do  this,  we  set  rg=l  and  we 
see  from(14)that  we  must  express  D,(c)-U,(c)  in 
terms  of  the  measured  signals.  The  :irst  layer  can 
be  depicted  as  in  Fig.  2.  Observe  tliat  (4)  becomes 

dj(t+T)  = -Ui(t)+2da(t).  (15) 

From  (11,  12,  13),  we  can  evaluate  the  difference* 
term  Dj(c)-U(e), 

Di(  e)-Ui(  e)  =4  P(  €)  = f'^*f[2d,(t)-Ui(t)] 

• [2do(t+€)-u^(t+€)]-uJt)u^(t+c)]dt  (16) 

or 

P(c)  = r_^%d,(t)d,(t+£)dt-  r_^j2d,(t)Uj(t+e)dt 


-|J'j2Ui(t)d,(t+ €)dt:  (17) 

hence,  P(e)  can  be  evaluated  from  a knowledge  of 
d,(t)  and  u^(t)  for  any  desired  s.  Observe,  also, 
that(14)with  i = l can  be  written  in  terms  of  P(£), 
using  ( 16),  as 

D,*i(€)  = P(e).  (18) 

We  should  point  out  at  this  stage  that  the  quan- 
tity u^(t)  needed  in(17)  is  only  available  through  the 
observation 


y(t)=Ui(t)+n(t)  (19) 

where  n(t)  is  the  additive  noise.  Consequently, P(e) 
is  not  physically  available:  however,  we  can  define 
P(e)  by  replacing  y(t)  for  Uj(t)  in  (17), 

P(€)  = j'_^“4d,(t)d,(t+e)dt  - 2d,(t)y(t+c)dt 

- r_^"2y(t)d,(t+e)dt  (20) 


Because  of  the  range  of  integration  in  (11),  we  can 
also  express  Dj(£)  as  Dj( c)  = J^^dpt+Tldjlt+T+cldt. 

We  use  this  form  of(ll)in  our  development  of 
D.je)-U(c). 


3 


which  ca.Q  also  be  writtea  as 

P(€)  = P(€)  + N(c)  , (21) 

where 


Note  that  in  general  the  a,  are  functions  of  d,^^(t),a 
signal  which  is  not  determinable;  however,  we  will 
show  in  the  following  section  that  when  the  input 
wavelet  is  narrow  eoough(not  necessarily  a spike), 
(28)  has  a unique  solution  for  the  reflection  coeffi- 
cients in  terms  of  observable  data. 


The  statistics  of  noise  term  N(c)  can  be  determined 
in  terms  of  those  of  n(t).  Using  P(e)  in(l8)  yields 

D,*i(e)  = IP  ?(£)  (23) 

whore  ^(e)  is  a known  quantity.  Eqaation(23)  is  a 
fundamental  relationship  which  will  be  used  in  the 
derivation  of  the  inverse  procedure. 

o.  derivation  of  the  normal  equations 

Equation  (7)  is  the  main  relation  which  will  be 
used  to  derive^the  reflection  coefficients.  Dividing 

both  sides  by  ^t,  and  identifying  the  resulting  coef- 
ficients by  So,  ....  5,.  we  get 

Sod,,j[t+KT]+3id,^i[t+(K-2)T]  r • • • 

9.-id,.i[t-(K-2)T]+e«d,,Jt-KT]  = do(t)  . (24) 

We  compute  the  coefficients  of  this  equation  by 
means  of  the  following  "least  squares"  criterion: 

a . T.l"  f 3od,  Jt+Kr]  + • ■ • + id.^Jt-Kr] 

-do(t)}»dt.  (25) 


The  result  of  this  minimization  is  equivalent  to 
multiplying  both  sides  of(7)  byd,*Jt+(K-2i)T]and  in- 
tegrating from  -•  to  • for  i = 0,  1, . . . , K.  By  either 
approach,  we  obtain  K+1  simultaneous  equations, 
which,  using  (11),  become 


I D,.,(0)  D,.„(2t)  ••• 

I Ai.i(2t)  D.,i(0)^-*- 
! : : ’ . 
IP».i12Kt)  • • • 


D,,,(2K:)]ri' 

1 

• iK 

I • 

• D.*t((J)  JW 


K 

2llt, 

0 


oil 

as 

a«*i 


(26) 

where  we  have  substituted  to=l+ro=2  to  represent 
the  marine  environment  and 


= j''^“do(t)d,..j(t+KT-2iT)dt,  i=0,  1,2 K. 

Substituting  for  D,,j(e)  from  (23),  we  find  that 
(26)  reduces  to 


7(0)  ■P(2t)  •••  P(2K-) 

1 

P(2t)  F(0)  ••• 

K 

OTa 

1 • • 

1 • • • • 

^(»T)  ' . . • Vlo) 

= 2nt,' 
0 

Note  that  the  (K+l)  x(K+l)  matrix  on  the  left  has  the 
Toeplitz  structure.  The  terms  N(0),  N(2t).  . . . 
which  appear  in  P(0),  P(2t),  . . . , are  random  vari- 
ables with  known  statistics;  they  will  be  zero  if  the 
seismic  data  is  noise  freed,  e. , o(t)eO).  Observe, 
also  that  the  first  and  last  elements  of  the  vector 
on  the  left-hand  side  of  (28)  are  unity  and  r,,  re- 
spectively, by  virtue  of  the  property  which  we 
stated  for  the  d,^^  causal  functional  equation. 

Equation  (28)  provides  the  second  point  in  our 
procedure  for  Identifying  the  reflection  coefficients. 


* 7.  SPECIAL  CASE  OF  NARROW  WAVELET 

Let  us  now  consider  the  case  where  d^(t)  does  not 
extend*  beyond  2t,  i.  e. , 

da(t)  = 0 t<  0,  t>2T.  (29) 

Since  the  time  of  arrival  at  the  K“  interface  is  Kt 
and  the  time  of  arrival  of  the  first  reflections  is 
(K+2)t. 

'l.*i(‘)  = 0 t<KT  (30a) 

c 

= 2nt5do(t-KT)  KT<t3(K+2)T  (30b) 

0 

= more  complicated  terms  t>(K+2)'f.  (30c) 
Prom  (27,  30a,  30b,  and  29)  we  see  that 

ai=  r^"d;,(t)d,^^(t+KT)dt  = 2nt.  r^'^d|(t)dt,  (31) 

0 

aad  that 


“'♦i'L^r‘^(‘'‘^«*if‘+^T-2iT]dt  = 0 i=l,...,K.  (32) 

Note  now  that(28)  will  have  precisely  K+l  unknowns, 
K of  them  in  the  vector  multiplying  the  Toeplitz  ma- 
trix and  one  on  the  right-hand  side,  oj. 

Finally,  Normal  Equation (28)  can  be  written  in  a 
compact  matrix  form,  as 

(33) 

where  P^  is  a (K+l)  y (K  + l)  Toeplitz  matrix  with  the 
first  row  being^  [P(0),  P(2t),  , . . , P(2Kt)];  a,  is  a 
K+l  column  vector  with  first  and  last  elements  1 
and  r„  respectively:  and,  c,=col(8,.  0, 0. . . . , 0)  and 

3,  = 2 (l-r*)  r^V(t)dt. 

1 . i t -0  » 

The  Normal  Equation  (33)  can  be  solved  for  This 
only  produces  one  of  the  K reflection  coefficients, 
namely  r,.  We  will  show,  in  the  following,thatinthe 
case  of  the  marine  environment,  nested  within  (33) 
are  a set  of  normal  equations,  the  solutions  of  which 
produce  each  one  of  the  reflection  coefficients.  The 
absence  of  this  useful  property  in  the  non-marine 
case  renders  the  procedure  of  this  paper  inapplic- 
able in  that  case. 

Let  us  now  hypothesize  a j-layer  system  (i.  e. , 
the  basement  layer  is  the  j**)  consisting  of  the  top 
j-layers  of  the  above  K-layer  system (K  s j). Clear- 
ly, from  (33),  we  have 


(34) 

If  this  condition  is  not  satisfied,  we  canalways  de- 
convolve the  data  to  achieve  this.  Since  the  re- 
quirement here  is  not  to  deconvolve  down  to  an 
impulse  function  (only  (29)  has  to  be  satisfied),  this 
results  in  a more  practical  solution. 

^For  a narrow  wavelet  and  s 5 2,  the  calculation  of 
P(e)  simplifies  since  (20)  reduces  to 


P(£)  = -2j*^  d5(t)y(t+c)dt . 
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where  will  agaia  have  1 and  r,  as  first  and  last 
elements.  We  shall  now  show  that,  in  the  case  of 
the  marine  environment,  P,  is  a (jtl)  x (j+1)  Toeplitz 
matrix  composed  of  the  t(m  left  corner  of  P,;  i.  e. , 
its  first  row  is  given  by  [P(0),P(2t) P(2j-)], 

For  the  moment  let  us  ignore  the  additive  noise 
term  in{20),  1-et  us  denote  by  ut(t)  the  response  of 
the  j-layer  systemli.  e. , the  term  u,lt)  in  Fig,  2 is 
replaced  by  u;(t)).  In(17),due  to  the ‘fact  that  d,(t)=0 
for  t>  2t,  the  last  value  of  u.(t)  contributing  to  P(e) 
is  u.(2T+e).  In  determination  of  P.,with  elements 
P(s),  £ = 0, . . . , 2jT,  for  a j-layered'svstem, therefore 
the  last  value  of  ul(t)  contributing  to  Pj  is  u;[2(j  + l)-]. 
On  the  other  hand,  u^(t)  is  the  response  of  the  K- 
layer  system,  and 

u,(t)  = uJ(t)  0sts2(j  + l)T  (35) 

since  the  first  return  from  the  interfaces  below  the 
j'*  will  not  a^ear  earlier  than  t=2(j  + l)T.  Hence. the 
elements  of  P,  which  are  functions  of  u-(t)  and  P, 
which  are  functions  of  u,(t),are  identica'l  when(35) 
is  satisfied.  In  other  words,  the  numerical  values 

oi  P(0) P(2j-)  will  be  identical  to  those  of  the 

K-layer  system  for  all  j s K.  Furthermore,  the  ad- 
ditive noise  term  in (21)  is  independent  of  the  num- 
ber of  layers,  as  is  evident  from  (22). 

The  set  of  normal  equations  given  by  (34),  for  j = 
1,  ....  K,  can  now  be  solved  for  the  vectors  a .,  and 
hence,  their  last  elements,  r,,  j=  I K.  The  ma- 

trix P^  is  Toeplitz  and  consequently,  the  Levinson 
algorithm  [l]  can  be  used  to  solve  for  the  vectors 
ij,  j = l,,..,K  recursively. 

Since  the  r.'s  are  reflection  coefficients,  for  the 
solution  to  this  problem  to  be  acceptable,  each  r 
must  be  less  than  unity  in  magnitude.  It  is  shown  in 
[9]  that  any  solution  of  (34)  with  I,  >0  for  all  j yields 
a set  of  r,'s  which  satisfy  this  condition.  More- 
over if  P«_is  positive  definite,  a compatible  solu- 
tion with  3,  > 0 is  guaranteed.  We  see  therefore, 
that  the  requirement  that  |r.|<l  has  nothing  to  do 
with  a specific  method  of  solution  of  the  normal 
equations  (i.e. , the  Levinson  procedure).  This  re- 
sult is  different  from  the  comparable  result  in[l- 
3],  where  one  is  left  with  the  impression  that  a 
specific  method  of  solution  leads  to  |rj|<l, 

8,  EXPERIMENTAL  RESULTS 

A seven  layer  system  was  chosen  with  the  fol- 
lowing reflection  coefficients;  r5=l:r.=0.1;  r,= 

0.15;  r3=-0,3:  r^=0.25:  rjsO.  12;  r,=0,‘05;  r,=0.  2. 

We  used  a one-way  layer  travel  time  and  data 
sampling  rate  of  20  msec  and  2 msec,  respectively. 

The  input  wavelet  was  chosen  as  depicted  in 
Figure  3.  Note  that  this  is  a non-minimum  phase 
function.  It  was  specifically  chosen  as  such  to  in- 
dicate that  the  method  introduced  in  this  paper  is 
not  limited  to  minimum  phase  wavelets.  Figure  4 
is  the  synthetic  seismogram  response  of  the  sys- 
tem. Figure  5 is  obtained  by  adding  white  Gauss- 
ian noise  with  variance  1 to  the  sampled  seismo- 
gram. Similar  results  were  obtained  for  vari- 
ances of  0.  1 and  10.  These  responses  were  then 
utilized  to  produce  estimates  of  reflection  coeffi- 
cient values.  The  results  of  a Monte-Carlo  simu- 
lation for  100  different  samples  of  noise  appear  in 
Tables  1 and  2.  Table  1 presents  the  mean  value 
of  the  estimates  and  Table  2 presents  the  estima- 


tion error  variance.  As  seen,  the  results  are  ex- 
act for  zero  noise  variance  (as  expected)  and  are 
quite  good  for  variances  of  0.1  and  1.  For  noise 
variance  of  10,  although  the  average  is  not  poor, the 
error  variance  Indicates  that  the  estimates  are  not 
very  reliable. 

A comparison  between  the  procedure  of  this  pa- 
Ppf  *tid  the  standard  procedures  described  in(l-4] 
is  warranted.  Let  us  consider  the  noise  term  n(t) 
to  be  white.  Clearly,  (22)  indicates  that  the  random 
variables  N(.)  have  finite  variances.  For  this  case 
I nit)  white],  had  we  performed  the  necessary  decon- 
volution and  sampling  required  by  the  classical  ap- 
proach to  the  inverse  problem,  the  resulting  N(>) 
random  variables  would  have  infinite  variance, clear- 
ly rendering  the  approach  meaningless.*  Of  course, 
"approximate"  deconvolution  will  eliminate  this 
problem  but  at  a great  sacrifice  in  the  information 
available  within  the  seismic  data.  It  should  also  be 
noted  that,  for  the  narrow  wavelets,  no  deconvolution 
is  required  by  the  procedure  outlined  in  this  paper, 
9.  CONCLUSIONS 

We  have  developed  a procedure  for  extracting 
reflection  coefficients  from  noisy  data  which  we 
feel  is  a substantiai  generalization  of  similar  pro- 
cedures which  have  been  reported  in  the  literature. 
Associated  with  these  earlier  procedures  are  Stan- 
dard Assumptions  and  Steps  (see  Introduction,  p.  1) 
which  include  requirements  that  the  data  be  noise 
free  and  that  the  observed  seismic  data  be  decon- 
volved. The  procedure  of  our  paper  avoids  these 
restrictive  requirements.  Furthermore,  our  proce- 
dure  totally  avoids  the  concepts  of  z-transforms, 
minimum  phase,  spectral  factorization,  forward  and 
reverse  polynomial  manipulations,  etc. , which  ap- 
pear in  the  literature  on  this  subject.  Finally.since 
our  derivation  is  so  straightforward,  it  suggests  a 
number  of  extensions,  including  the  following,  which 
are  presently  under  study:  (1)  nonstandard  locations 
of  source  and  sensors  (e.  g. , both  in  the  first  layer): 
(2)  minimum  mean-square  estimation  in  the  non- 
marine environment:  and  (3)  optimal  prefiltering  of 
noisy  data. 
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